Flow stabilization with active hydro dynamic cloaks 
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We demonstrate that fluid flow cloaking solutions based on active hydrodynamic metamaterials 
exist for two-dimensional flows past a cylinder in a wide range of Reynolds numbers, up to approxi- 
mately 200. Within the framework of the classical Brinkman equation for homogenized porous flow, 
we demonstrate using two different methods that such cloaked flows can be dynamically stable for 
Re in the range 5-119. The first, highly efficient, method is based on a linearization of the Brinkman- 
Navier- Stokes equation and finding the eigenfrequencies of the least stable eigen-perturbations; the 
second method is a direct, numerical integration in the time domain. We show that, by suppressing 
the Karman vortex street in the weekly turbulent wake, porous flow cloaks can raise the critical 
Reynolds number up to about 120, or five times greater than for a bare, uncloaked cylinder. 
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Metamaterials (MMs), or artificially structured com- 
posites, have been proposed as a supplement to natural, 
molecular media that extends material properties into 
the domains not covered by naturally formed substances. 
Supplemented with a macroscopic design methodology, 
known as transformation optics (acoustics, and so on), 
MMs enable a new wave applications by offering extra 
flexibility in manipulating the dynamics of waves and 
matter. Recently, it was proposed [1 that hydrodynamic 
metamaterials (HDMM) consisting of a fluid-filled solid 
matrix with a variable fluid permeability can enable ex- 
otic flow regimes, in which pressure and velocity differ- 
entials are confined to a small volume occupied by the 
permeable shell. These flows were demonstrated in three 
dimensions, for a spherical object surrounded by a con- 
centric porous shell with anisotropic, and partially nega- 
tive, permeability. As explained below, such a property 
would necessitate the use of active hydrodynamic meta- 
materials, defined as HDMMs with active elements that 
can accelerate the fluid; the latter elements can be de- 
vised, for example, from electrically powered and con- 
trolled micro-pumps [2j [3] . 

The wider concept of active, or externally powered, 
MMs is currently a subject of fundamental studies in 
the fields of electrical engineering, optics, and acous- 
tics. Here, before attempting to propose a physical im- 
plementation of an active HDMM, we wish to address the 
question whether such media would be of any benefit in 
hydro- or aerodynamical engineering. This paper serves 
to demonstrate that active HDMMs are capable of lever- 
aging the critical Reynolds number (Re) at which the 
flow past an object transitions from laminar to weakly 
turbulent. 

In comparison with the three-dimensional flow past a 
sphere, the 2D case is substantially more challenging. 
From the theoretical perspective, exact analytical solu- 
tions for flow past a cylinder are not possible even in the 



limit of small Re; Oseen's equation can be used to obtain 
the lowest-order approximation [4 . Thus, 2D flows are 
more sensitive to the approximate treatment (or neglect) 
of the momentum advection term ((u-V)u) in the Navier- 
Stokes equations. This term is responsible for chaotic 
motions that are bound to grow from small perturbations 
when Re exceeds a certain critical value; well-developed 
stochasticity of a flow with stationary boundary condi- 
tions is known as turbulence g]. 

For 2D flows past a cylinder, the growth of perturba- 
tions is a concern even at low Re (less than a hundred), 
where the spontaneous formation and separation of os- 
cillating eddies in the downstream wake gives rise to the 
Karman vortex street [5 . Vortex shedding has lots of un- 
desirable practical consequences, including vibrations of 
long fluid-protruding objects (towers, chimneys, aerials) 
and bubble generation in gas-saturated liquids (subma- 
rine periscopes). Engineering effort was made to control 
the wake and suppress vortex shedding by manipulating 
the separation of vortices from the boundary layer by 
both passive means such as fins and other vortex spoil- 
ers [Bj, as well as active devices [5j[7J[E]- In this paper, 
we propose that fluid-permeable active flow "cloaks" can 
suppress the formation of eddies in a wide range of Re 
numbers, thus leveraging the onset of turbulence. 

We begin our analysis by postulating the time- 
dependent nonlinear Brinkman equations for the macro- 
scopic description of flows [9HT2] in HDMMs: 
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We wish to simplify the analysis by setting the poros- 
ity parameter e, which represents the filling fraction of 
the fluid, to a constant value of unity. This approxima- 
tion is consistent with the expectation that the desirable 
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FIG. 1. (color online), (a) Schematic of the structure and 
the problem setup, (b) Permeability profiles designed for 
flows with various flow rates, labeled by the corresponding 
Reynolds number Re = puoa/p. Steady-state cloaking solu- 
tions: distribution of the vertical velocity u y normalized to 
uo (shade) and streamlines of flow (black lines) for solutions 
optimized at various flow rates: (c) Re = 0.5, (d) Re = 190. 

porous media will have a very small filling fraction of 
the solid phase, which would allow them to have a high 
permeability (low resistance to the flow). Generalization 
of the forthcoming analysis to specific values of e < 1 is 
straightforward. 

With this simplification, the steady-state form of the 
Brinkman equation for an incompressible fluid (p = 
const) with a constant free-stream dynamic viscosity 
p and coordinate-dependent inverse permeability tensor 
ft -1 becomes 

p(u ■ V)u = — Vp + pS7 2 u — /ift _1/ a, (3) 

combined with the steady-state continuity equation, 

V ■ ?1 = 0. (4) 

Within the framework of Eq. ([2]), we solve the problem 
of two-dimensional flow cloaking [1 , or wake elimination. 
We then proceed to demonstrating that wake elimina- 
tion indeed results in flow stabilization at flow rates that 
normally result in turbulence. The schematic of the ge- 
ometry under consideration is shown in Fig. [TJa). The 
boundary conditions are no-slip (u = 0) at the surface of 
an impermeable cylinder (r = a), and plug flow in the 
y-direction (u = u^y) at infinity (r — >• oo). Elsewhere, in- 
cluding the transition from porous to free flow domains 
at r = 6, the continuity conditions [13] are assumed for 
the velocity field and for the fluid stress tensor, includ- 
ing its pressure and shear components. Whenever a spe- 
cific calculation occurs, we assume a = 0.5 mm, b = 2a, 



p = 1 • 10 -3 Pas, and p = 1 g/cm 3 . When referring 
to normalized permeability, the normalization constant 
fto = a 2 /4 is assumed. 

Previously, approximate solutions were reported for 
the case of flow past a spherical cloak pQ. The proposed 
technique was based on the introduction of the (Stokes) 
stream function and solving the linear equation in the 
regime with vanishingly small Re. However, for the flow 
past a cylinder, it is well known that the linear Stokes 
equation cannot give meaningful solutions that satisfy 
the boundary conditions at r = a and infinity; the origin 
of the problem is the existence of a logarithmic funda- 
mental solution. This issue can be partially alleviated 
by considering Oseen's equation, which uses an approx- 
imate, linearized form of the advection term instead of 
neglecting it entirely. Still, this term, due to its angular 
symmetry, would break down the separation of variables 
in the polar coordinates, thus rendering the method of 
Ref. pQ much more difficult to apply. Thus, we choose to 
obtain the wake-free (cloaking) solutions numerically by 
employing a gradient-based optimization technique. 

It was shown in Ref. [1 that wake-free solutions ex- 
ist when the permeability is spherically symmetric but 
locally anisotropic with the spherically-uniaxial type of 
anisotropy. Here, we show that approximately wake- free 
designs can be obtained by assuming locally isotropic 
and cylindrically-symmetric permeability n(r). For solv- 
ing the forward steady-state problem with any given 
ft(r), our approach employs a standard two-dimensional 
Galerkin (finite element) method provided by simulation 
software COMSOL Multiphysics [14], specifically, in its 
"free and porous flow" interface, which includes both 
Brinkman and Navier- Stokes equations. The simulation 
domain throughout this work is a rectangle of width 
W = 20a horizontally and height H = 60a vertically. 
The inflow boundary is positioned at y = —H/3 and the 
outlet is at y = 2i^/3, with the center of the cylinder 
at the origin. The forward problem is then embedded 
into a gradient-based optimization algorithm (SNOPT), 
which is also part of the same program. This integrated 
forward- inverse solving strategy lets us use the built- 
in semi-analytic sensitivity analysis, which produces the 
gradient of the optimization goal with respect to all opti- 
mization variables in one step, which is then used by the 
gradient-assisted optimization routine. The optimization 
goal is to minimize the norm of velocity deviation from 
a perfectly uniform plug flow, i.e. 

W = J \u-u \ 2 dV, (5) 

where the integral is taken over the exterior of the cylin- 
der of radius r — b which represents the "cloak". The 
quantity W can be seen as a measure of wake in the 
flow past the structure. The optimization variables are 
the nodal values of the unknown permeability function, 
ft(r), discretized on a certain finite element mesh. These 
nodal values are then extruded and interpolated onto the 
entire annulus a < r < b. 



In this fashion, we obtain a fully converged solution 
with the flow rate corresponding to Re = 0.5, and then 
use that solution as an initial guess for the next run of 
the optimization solver at a slightly higher inflow velocity 
uq. This procedure yields fully converged (in the sense of 
the gradient-based optimization) solutions, showing very 
little velocity perturbation ([5| in the exterior of the per- 
meable cylinder. Several such solutions are illustrated in 
Fig.Qb-d). 

A key feature of the obtained solutions is that they 
contain regions where n~ x < 0. While we are not aware of 
physical implementations of negative n as passive porous 
media, we may conceive it to be implemented as an active 
HDMM. The latter is defined here as a fluid-permeable 
structure that applies an accelerating force to the fluid. 
To understand how an active HDMM can implement the 
proposed solutions, one should realize that the Brinkman 
Eq. ([3]) is nothing more than the momentum transport 
equation with three volumetric forces: pressure gradient 
force, viscous shear force and, lastly, the force exerted by 
the solid components of the HDMM upon the fluid, 

F s = —fin^u. (6) 

In passive media (or media with n > 0) F s is typically 
contra-directed with the average flow direction u, so that 
it does work against the flow (F s • u < 0) and thus de- 
celerates it. However, in an HDMM an embedded array 
of thrust-generating elements (such as mini- or micro- 
pumps [21 [3]) may generate a force F s with an arbitrary 
direction relative to {?, including F s • u > 0. Therefore, 
our HDMM cloaking solutions are to be implemented in 
practice as follows. Once the spatial distributions n(f) 
and macroscopic fluid velocity u(f) are found using the 
technique presented here, negative permeability regions 
are replaced with an equivalent body force, f thrust, us- 
ing the trivial relationship pn~ l (r)u(f) = pn~ x (r)u(f) — 

fthrustif), where n p (r) > is an arbitrary positive func- 
tion of radius. We emphasize that the origin of negative 
permeability in our results is the ansatz of Eq. ^ used 
to express the volumetric force needed to sustain the de- 
sired flow pattern (F s ) in terms of the macroscopic fluid 
velocity (u). This ansatz is used merely for historical 
reasons, and it is neither expected to be valid at high Re 
nor required to obtain the types of flows reported here. 

In order to address the question of dynamic stabil- 
ity of the stationary solutions obtained, we first employ 
the method of small perturbation [15]. The technique 
is based on the linearization of the full time-dependent 
Eq. Q with respect to the steady-state solution (u 3 , p s ) 
of the nonlinear Eq. Assuming that the velocity and 
pressure perturbations (v, p') depend on time as e _A£ , 
the following eigenvalue problem can be formulated: 

p(u s -V)v + p(v-V)u s + Vp' — pV 2 v + pn~ x v — \pv, (7) 

which must be supplemented by the linearized continuity 
equation, V • v = 0. 
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FIG. 2. (color online). Stability analysis using harmonic 
perturbation method for several fluid cloaks. Distribution 
of velocity perturbation (shade) and its streamlines (black 
lines) for the eigenmode with the lowest Re(X); cloaks with 
(a) Re = 1 (unstable), (b) Re = 5 (stable), (c) Re = 100 
(stable), (d) Re = 120 (unstable). 



Re 


Ai 


A 2 


A 3 


0.5 


-4.093 


0.3456 ± 0.05762 


0.4466 ±0.13932 


1 


-2.621 


0.6843 ± 0.19142 


0.818 ± 0.46842 


4 


-0.4757 


3.008 ± 0.70272 


3.4721 ± 2.4045i 


5 


0.3651 


3.823 ± 0.59222 


4.477 ± 3.0082 


10 


3.369 


7.4493 


8.7070 


50 


18.618 ±40.7152 


48.461 ± 55.73H 


51.655 


100 


9.748 ± 125.922 


93.835 


94.886 ± 96.4842 


119 


8.950 ± 154.3872 


113.795 ± 114.2672 


114.164 


120 


-11.718 ± 134.1772 


115.740 


118.602 ± 73.0882 


190 


-21.490 ±211.9382 


85.055 


182.137 



TABLE I. The three lowest eigenvalues (ordered by increasing 
Re A) for cloaks designed for different Re. 



Once the eigenvalue problem is solved, the eigenvalues 
A n can be used to determine whether the solution is sta- 
ble with respect to small perturbations. The steady-state 
solution is unconditionally stable when all eigenvalues A n 
lie in the right half-plane, i.e. Re A n > Vn. Thus, the 
critical velocity at which the transition from laminar to 
turbulent flow occurs can be determined by looking at 
A = min n ReA n as a function of 220; once A becomes 
negative, the flow is no longer stable. 

First, we calibrated this method by solving the clas- 
sical problem of laminar flow past a bare, impermeable 
cylinder of radius a (results not shown). The spectra we 
obtain indicate that such flow is unconditionally stable 
for Reynolds numbers Re = puoa/p up to Re^r = 23.5, 
and bears unstable eigenmode (s) above that Reynolds 
number. The instability-causing eigenmode in the latter 
regime has the appearance of the familiar Karman vortex 
street. This finding is in general agreement with other 
theoretical and experimental results [I6j [17] . 

We then apply this method to porous cloaks with var- 
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ious flow rates, which are labeled by the correspond- 
ing Reynolds number of the flow, Re = puoa/fi, cal- 
culated with respect to the "cloaked" cylinder radius 
a. The lowest three eigenvalues for cloaks with vari- 
ous Re are listed in Table |TJ Noticeably, an unstable 
eigenmode exists in the range Re < 4, which becomes 
stable at Re « 5 and remains stable at higher Re. A 
snapshot of this eigenmode's velocity perturbation pro- 
file is shown in Fig.[2|a,b), both below (a) and above (b) 
the transition Re. This behavior — flow stabilization 
with increasing flow rate — is drastically different from 
the usual laminar- turbulent flow transitions. Apparently, 
this form of instability can be attributed to the presence 
of an active (negative-permeability) medium, which adds 
momentum to the flow in proportion to the local velocity 
(see Eq. (|6| and the discussion below it). Not surpris- 
ingly, a hypothetical medium with a momentum pump 
whose strength grows with the flow rate acts as an am- 
plifier of local fluctuations, causing the system to exhibit 
an instability even at very small Re well within the Stokes 
flow regime. 

In the range Re = 5 - 119 (see Tableland Fig.|2fc)), 
all eigenmodes existing in our simulation domain are sta- 
ble, which leads us to believe that the steady-state solu- 
tions obtained in this regime are unconditionally stable. 
This conclusion is further confirmed below by another 
method. At Re = 120, an eigenmode corresponding to 
the formation of an oscillating Karman vortex street past 
the cylinder (Fig. [2^d)) becomes unstable, and remains 
unstable for all higher Re that were simulated (up to 
Re — 192). One may observe from Fig. [2^d) that the 
eigenmode leading to flow instability in that regime ap- 
pears to be seeded by the perturbation in the small region 
past the cylinder near its back surface, where the velocity 
deviation (or wake) is not completely eliminated in the 
stationary solution (see Fig. []Jd)). Since the cloaking 
structures obtained through numerical optimization of 
the permeability deviate noticeably from a perfect wake- 
free goal ([5| at Re > 100, one can therefore expect that 
a more accurate cloak with a much smaller value of ([5| 
would also possess a more stable perturbation spectrum. 

To validate the findings of our perturbative stability 
analysis, we perform full transient simulations based on 
the full time-dependent Eq. (2|. The initial value prob- 
lem was solved for cloaks with several u$, with the same 
boundary conditions as were used for the stationary anal- 
ysis, with the initial values u(x,y,t = 0) = u^y and 
p(z,y,t = 0) = 0. 

For the cloak designed to operate at Re = 1, which, 
for the cylinder dimension 2a = 1 mm corresponds to in- 
flow velocity uq = 2 mm/s, we observe that the solution 
quickly (within 0.1-0.2 s) reaches the steady-state solu- 
tion that was found in the stationary analysis (Fig. [3^a- 
b)). However, as evident from Fig. |3ja), at t ~ 5 s since 
the simulation start, the velocity in the vicinity of the 
structure begins to grow exponentially, and at t = 6.6 s 
(Fig.|3jc)) the velocity profile matches precisely with the 
unstable eigenmode predicted by the perturbation analy- 
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FIG. 3. (color online). Time-dependent analysis for fluid 
cloaks optimized at different Re. Cloak operating at Re — 1: 
(a) vertical velocity on the surface of the cloak at position 
x = 6, y = vs time (log scale); (b) v y velocity component 
profile at t = 0.3 s (before instability) and (c) t = 6.6 s (after 
unstable eigenmode growth). At Re = 100: (d) vertical veloc- 
ity at a surface point vs time; (e) v y snapshot at t = 0.0015 s 
(before stabilization) and (f) t = 0.15 s (after stabilization). 



sis (see Fig.|2ja)). Pure exponential growth of this eigen- 
mode continues throughout the end of the simulation pe- 
riod (10 s). 

Similar behavior, but on a longer temporal scale is ob- 
served for the cloak designed to operate at Re = 4 (not 
shown). As seen from Table [IJ the unstable eigenmode is 
characterized by the negative real part Re A = —0.4757, 
which is less by factor 0.18 than Re A = —2.62 of this 
eigenmode at Re = 1. Thus, one may expect that 
the time needed for the instability to develop is about 
1/0.18 « 5 times greater than in the previous case. This 
prediction is confirmed by the transient simulation: the 
instability begins to grow at t ~ 25 s (not shown). 

The cloak designed for Re = 100 (inflow velocity 
uq = 0.2 m/s) shows a drastically different behavior. The 
flow reaches steady state in about 0.04 s since initializa- 
tion time, and remains completely unchanged through- 
out the rest of the simulation period T = H/uo = 0.15 s, 
which is chosen such that the fluid starting at the in- 
let has enough time to travel throughout the very long 
simulation domain. This confirms that the steady-state 
solution found above is dynamically stable. Virtually 
identical behavior is observed in the flow with Re = 119 
(not shown), which is the upper limit of the dynamically 
stable regime in our calculations. 

Above the transition (Re = 120), numerical integra- 
tion of time-dependent equations reveals an instability 
(Fig. [I]). The spatial profile of velocity disturbance, and 
its oscillating nature allows one to immediately identify 
the vortex sheet eigenmode shown in Fig. [5|d). This 
behavior is seen in the entire range of velocities studied 
(Re = 120 — 192). The growth rate of this instability dur- 
ing its unsaturated (exponential) growth, as well as its 
oscillation frequency are consistent with correspondingly 
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FIG. 4. (color online). Transition to turbulence in the de- 
signed cloaks. At Re =190: (a) vertical velocity at a surface 
point (x = b, y = 0) vs time; (b) v y profile at t = 0.004 s (be- 
fore instability growth) and (b) t = 0.4 s (after considerable 
growth of instability). 



the real and the imaginary parts of the lowest eigenvalue 
(Ai) shown in Table [TJ 

The close agreement between the two independent sta- 
bility analysis methods is evidence, that, within the ap- 
plicability domain of the Brinkman-Navier-Stokes equa- 
tions, the fluid-permeable structure described above is 
able to increase the critical velocity u cr , which marks 
the transition from laminar to turbulent flow, by factor 
of five. According to the law of hydrodynamic similar- 
ity, u cr is a function of the object diameter, 2R (which 
also equals its cross-section in two dimensions), kinematic 
viscosity v = /i/p, and a single dimensionless coefficient, 
Re cr which itself is a function of only the object shape: 



u cr = He ^ v . Thus, the increase in u cr can be stated ei- 
ther as an increase in Re cr , the critical Reynolds number, 
or as a reduction of effective diameter, 2R e R , defined ac- 
cording to u cr = R ^ - , where Reir = 23.5 is the critical 
Re of the bare cylinder. In that sense, the fluid sees the 
cylinder of radius R surrounded by the cloaking structure 
as a cylinder of radius i? e ff, which can be substantially 
smaller than R, according to the findings above. 

In conclusion, we have demonstrated that approximate 
fluid flow cloaking solutions exist for flows past a cylin- 
drical obstacle. These solutions require a mixture of pos- 
itive and negative permeability inside a cylindrical shell, 
which can be locally isotropic. We show that negative 
permeability gives rise to dynamic instability for flows 
with Reynolds numbers Re < 5. However, a family of 
solutions exist for all Re in the range 5 — 119 where the 
steady-state flow is unconditionally stable. This find- 
ing hints at the possibility of maintaining laminar flows 
around cylinders at velocities five times the critical veloc- 
ity for a bare cylinder. Further minimization of the wake 
at high Reynolds numbers should lead to dynamically 
stable solutions at the correspondingly higher velocities. 
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